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ABSTRACT 


A numerical  approach  to  solving  the  two-dimensional , 
incompressible  Navier-Stokes  equations  is  presented 
and  is  used  to  study  the  stability  and  evolution  of 
disturbances  in  a boundary  layer.  Previous  numerical 
approaches  have  used  streamfunction-vorticity  formu- 
lations that  take  advantage  of  the  special  properties 
of  two-dimensional  flow.  The  present  method  solves  a 
finite-difference  form  of  the  momentum  equations  and 
may  therefore  be  extended  to  three-dimensions  more 
readily.  Computations  using  this  technique  are 
compared  with  the  result  of  linear  stability  theory 
for  small  amplitude  disturbances  and  are  found  to 
give  satisfactory  agreement,  while  calculations  at 
large  amplitude  support  the  conjecture  that  non- 
linear effects  can  be  destabilizing. 
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INTRODUCTION 


One  of  the  principal  interactions  between  a fluid  and  a solid 
structure  is  the  tangental  force  due  to  relative  motion.  At 
low  velocities,  when  the  flow  is  still  laminar,  the  drag  force 
is  small,  but  at  higher  speeds,  when  the  flow  has  become 
turbulent,  the  drag,  as  well  as  the  heat  transfer,  increases 
enormously.  It  is  of  practical  interest  to  be  able  to  under- 
stand the  mechanisms  leading  to  the  transition  from  a laminar 
to  turbulent  flow  and  to  be  able  to  predict  and  control  that 
phenomena. 

The  onset  of  turbulence  is  preceeded  by  a motion  in  the  boundary 
layer  that  is  certainly  three-dimensional  and  probably  highly 
nonlinear.  Any  theory  that  would  explain  the  transition  process 
must  address  these  two  issues  in  some  manner.  At  its  initial 
stages,  however,  transition  appears  to  involve  fairly  regular, 
small  amplitude,  two-dimensional  disturbances  that  are  capable 
of  extracting  energy  from  the  mean  flow  and  thereby  grow 
rapidly.  This  motion  is  within  the  power  of  analysis  to 
investigate  as  an  eigenvalue  problem  and,  consequently,  it  has 
been  extensively  studied. 

The  extension  of  these  results  to  true  initial  value  problems, 
large  amplitudes,  and  three- dimensions  has  proven  to  be  quite 
difficult  for  analysis  and  so  it  is  desirable  to  develop  other 
approaches  to  assist  in  understanding  these  effects.  One 
promising  approach  is  the  direct,  numerical  solution  of  the 
governing  equations.  Encouraging  results  have  already  been 
obtained  by  Easel  (1976)  and  Murdock  (1977)  using  a stream- 
function-vorticity  formulation  of  the  equations  of  motion. 

Easel  employed  a three-time-level  finite-difference  solution 


algorithm,  while  Murdock  has  combined  an  orthogonal  function 
expansion  in  one  spatial  direction  with  a finite-difference 
formulation  in  the  other.  The  two  methods  appear  to  give  good 
agreement  with  linear  analysis  and  with  one  another. 

The  present  work  explores  another  numerical  approach;  the  finite- 
difference  solution  of  the  momentum  form  of  the  Navier-Stokes 
equations.  Three  variables  (two  velocities  and  a pressure)  must 
now  be  calculated  instead  of  two  (a  streamfunction  and  vorticity), 
however,  the  numerical  technique  is  readily  extended  to  three- 
dimensions  where  the  concept  of  a streamfunction  is  no  longer 
valid.  The  finite-difference  algorithm  used  differs  from  that 
of  Fasel  in  two  respects.  The  present  scheme  introduces  no  numer- 
ical viscosity  into  the  solution,  while  Fasel 's  scheme  is 
slightly  dissipative.  Secondly,  the  present  scheme  requires 
storing  the  flow  information  at  only  one  time  level,  instead 
of  at  two  levels  as  in  Fasel ' s approach.  The  consideration  of 
computer  storage  may  be  a critical  one  when  three-dimensional 
calculations  are  attempted.  Both  schemes  involve  systems  of 
nonlinear  algebraic  equations  which  must  be  solved  Iteratively. 

The  work  presented  here  describes  the  details  of  the  problem 
formulation  and  numerical  solution  procedure  as  applied  to 
spatially  periodic,  parallel  flows.  Comparisons  with  the  results 
of  linear  analysis  are  made  and  verify  the  correctness  of  the 
numerical  solutions.  Finally,  the  effects  of  large  amplitude  on 
subcritical  stability  are  investigated  and  indicate  that  non- 
linear effects  can  be  destabilizing. 


2.  THE  MATHEMATICAL  FORMULATION 


The  motion  of  the  fluid  is  assumed  to  be  described  by  the  two- 
dimensional,  incompressible,  time-dependent  Navier-Stokes 
equations.  Normalizing  the  spatial  distances  by  the  momentum 
thickness  of  the  boundary  layer,  6,  the  velocities  by  the  free 
stream  value,  U^,  and  time  by  6/U^,  the  equations  of  motion 
become 


3u  . , 3u  . 9u  . 3p  ln2 

5T  * “3^  + ''W  ^ r’  “ 


3v  . 9v  , 3v  . 3p 
3t  3x  3y  3y  R 


+ 1^  = 0 

3x  3y 


Here,  x and  u are  the  coordinate  and  velocity  in  the  streamwise 
direction,  y and  v are  the  coordinate  and  velocity  in  the 
direction  normal  to  the  boundary,  p is  the  fluid  pressure  divided 
by  the  mass  density,  and  R = U^6/v  is  the  Reynolds  number,  where 
V is  the  kinematic  viscosity.  The  Laplacian  operator  is 


V2  = 


It  is  convenient  to  define  a base  flow,  denoted  by  capital 
letters,  that  satisfies  the  steady  flow  equations 


u|U  + v|U  + |P  = iv^u 

3x  3y  3x  R 


3 


^9x  ''ay  3y  R 


. 8V  ^ 
3x  ay 


For  high  Reynolds  number  flows,  these  equations  can  be  solved 
either  directly  or  by  using  boundary-layer/potential  flow 
techniques.  Disturbances  to  the  base  flow  then  satisfy 


9v:  + uii:  ^ .3v  ^ .3v:  ^ ^ ^ 

at  ^ 3x  ^ 3x  ^ 3x  ay  3y  ay 


+ Ip:  = 

3y  R 


+ 1x1  = 0 

TF  ^ ay  ^ ' 


The  boundary  conditions  appropriate  to  these  equations  are  that 
the  disturbance  quantities  vanish  far  from  the  solid  boundary, 

u',v',p'  -^0  as  y , (10) 

the  disturbance  velocities  vanish  on  the  solid  boundary, 

u^  = v'  = 0 at  y = 0 , (11) 

and,  by  evaluating  (8)  on  the  boundary  and  using  (9),  that 
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rnsmemm 


mmmm 


^ = _i  iV  at 

3y  R 3x3y 

Boundary  conditions  in  the  streamwise 
nature  of  the  problem  being  solved. 


y = 0 . (12) 

direction  depend  on  the 


The  techniques  developed  in  this  work  are  capable  of  solving 
the  disturbance  equations  for  the  most  general  base  flows,  but 
here  we  shall  consider  the  special  case  of  parallel  flow,  when 
U = U(y)  and  V = 0.  Strictly  speaking,  the  only  such  flow 
satisfying  this  condition  is  the  parabolic  distribution.  In  the 
thin  boundary  layers  that  form  at  high  Reynolds  numbers,  however, 
the  normal  velocity  field  is  extremely  small  and  the  assumption 
of  parallel  flow  is  nearly  met  (Schlichting,  1960,  pll6ff).  The 
disturbance  equations  can  then  be  simplified  to  the  form 


3u" 

3t 

, .dU  ^ 

+ V -T-  + 

dy 

^.3u' 

9y 

+ ^ .. 
3x 

= 

R 

3v' 

3t 

+ 

3y 

^ Y ■■ 

3y 

= 

R 

(13) 


(14) 


= 0 

3x  3y  ^ 


(15) 


The  linearized  form  of  equations  (13)  - (15)  have  been  extensively 
studied,  especially  the  growth  rates  and  spectral  properties  of 
periodic  disturbances.  The  assumptions  of  linearization  and 
periodicity  allows  the  equations  above  to  be  transformed  into 
a single  ordinary  differential  equation,  the  Orr-Sommerfeld 
equation.  The  resulting  body  of  knowledge,  known  as  linear 
stability  theory  (LST),  offers  a quantitative  test  of  the  direct 
numerical  solution  of  the  equations. 

The  results  discussed  in  this  work  are  for  parallel  flows  in 
which  U is  taken  as  the  horizontal  component  of  the  velocity 


over  a flat  plate,  the  Blasius  profile.  In  order  to  compare  with 
LST,  the  disturbance  flow  is  assumed  to  be  periodic  with  wave- 
length X in  the  streamwise  direction.  There  are  then  two  param- 
eters describing  small  amplitude  motion,  the  Reynolds  number,  R, 
and  the  dimensionless  wavenumber  a = 2it6/X. 
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3.  THE  NUMERICAL  SOLUTION  TECHNIQUE 


The  equations  (13),  (14),  and  (15)  subject  to  boundary  conditions 
(10),  (11),  and  (12)  are  to  be  solved  in  a region  0 £ x ^ X/6 
and  0 _<  y < ».  The  region  is  first  divided  into  finite  sized 
cells,  of  constant  size  in  the  horizontal  and  variably  spaced 
in  the  vertical  (see  Figure  1).  Within  each  cell,  the  horizontal 
velocity  is  defined  at  the  center  of  the  left  face,  the  vertical 
velocity  at  the  center  of  the  bottom  face,  and  the  pressure  at 
the  cell  center  (see  Figure  2) . 

The  following  finite-difference  form  was  used  to  represent  the 
time  rate  of  change  of  the  discrete  velocities 


(— ) 

^3t^i,j 


U.  (u,^,  , - u.  , .) 


26x  "j  ^“i+l,j  ^i-l,j 


46 


6x  [^i+l,j  ^'^i+l,j  ^i+l,j^ 

“1-1. d <“i,j  * “i-i.j>] 

4^  [“i.J+1  ^''i.d+l  * ''l-l.j+l' 

“l.j-1  (''l.d  * ''1-1, J>] 

1/dU. 

4^Hy^j  ^^^4,3  '^i.j+l  '^i-lJ  '^i-l,j+l) 


6x  ^Pi,j  “ Pi-1, j^ 


(continued) 


-jC  T ' ' ■ t - r-  , r r' 


Here  6y.  , = i(5y.  + 6y.  Also,  we  have  dropped  the  prime 

from  the  lower  case  quantities  since  no  confusion  with  total 
velocities  or  pressures  can  result.  The  difference  form  used 
for  the  nonlinear  terms  is  that  of  Piacsek  and  Williams  (1970). 

If  I and  J represent  the  total  number  of  columns  and  rows  used, 
with  columns  i=l  and  I and  rows  j=l  and  J being  boundary  cells 
outside  the  physical  domain  , then  the  boundary  conditions 
become 


u 


""i,J-l 


Pi,J  "Pi,J-l 


(18) 


''1,2  - ° • “1,1  = -“1,2  , 


J*!,!  Pi, 2 * R6x  ^“i+1.2  ' “i,2* 


u 


I.j 


= u 


2,j 


Pl.j  = P2.J 


(19) 


(20) 


u 


I.j 


= u 


I-l,j 


'l.j 


= V 


I-l,j 


Pi 


, J 


= p 


I'l.j 


. (21) 


The  pressure  field  is  determined  so  as  to  guarantee  that 
equation  (15)  is  satisfied.  Before  discussing  the  details  of 
how  the  pressure  is  calculated,  it  is  convenient  to  define  the 
tilda  finite-difference  operators 


. A/  ^ 

/ 9 U \ _ / 9 11  \ 1 f T\  r\  ^ 

^3ta,j  “ ^3t^i,j  " 6x  ^Pi,j  " Pi-1, j^ 


(22) 


/lY)  = ,iv. 


(p.  ^ 4_1 ) . (23) 


3t''i,j  ''3t''i,j  i.j  i.j-1 
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that  is,  they  are  the  time  rate  of  change  operators  without  the 
pressure  terms.  The  solutions  are  advanced  in  time  by  the  time- 
average,  implicit  algorithm 


where  the  superscript  indicates  the  time  level  and  6t  is  the  finite 
time  step.  This  scheme  when  used  with  central,  space  differencing 
introduces  no  numerical  dissipation  (Roach,  1976,  p.84).  The 
velocity  fields  are  then  given  by  the  formulae 

, n+1  _ n 1 r^Suxn  .^u.  T 6t  ,r„n+l  ^n+1 


_ % 5t  / n+1  n+1  . 

= “i,J  - 263  <Pi,J  - 


l.J 


6t  /_n+l 

26yj.j  <Pl,J 


n+1  . 


_ 6t  , n+1  „n+l  . 

'^i,j  " 26yj_^  ^^i.j  ” ^i,j-l^ 


Because  the  time  derivatives  involve  products  of  the  fluid 
velocities,  (26)  and  (28)  are  coupled,  nonlinear  algebraic 
equations.  Furthermore,  the  pressure  field  at  the  new  time 
level  must  be  chosen  so  that  the  new  velocities  are  divergence 
free,  that  is 


1 . n+1  n+1.  ^ 1 , n+1  n+1.  . 

<“1+1, j - "i,j>  ^ JTi  ‘"i,!!!  - ''i.j’  ' ° • 


This  last  relation  leads  to  the  following  linear,  elliptic 
boundary  value  problem  for  the  pressures 


1 r n+l  n^n+l  . n+1  1 

(5x)^  ~ 


1 r 1 , n+l  n+l. 

y,.  Uyj+i 


Sy 


1 /_n+l  „n+l  ^1 


_ 2 f 1 ,'v.  a, 

6t  L<Sx  ^'^i+l.j  ~ ^i,j^ 


1 

5yj  ^’^i,j+l 


- V.  .)]  (30) 

J.  , J 


The  solution  is  advanced  to  the  n+l  st  time  step  through  the 
following  iterative  scheme 


1.  The  new  velocities  and  pressures  are  set  equal 
to  the  previous  values. 

2.  The  tilda  velocities  are  computed  using  (25)  - 
(28). 

3.  The  system  of  equations  (30)  is  inverted  to 
obtain  a better  guess  of  the  pressure. 

4.  Equations  (26)  and  (28)  are  used  to  obtain 
a better  guess  of  the  velocities 

5.  Steps  (2),  (3),  and  (4)  are  repeated  until  the 
quantities  do  not  change  by  some  predetermined 
amount . 

The  solution  of  (30)  is  obtained  by  decomposing  and  backsolving 
the  banded  pressure  matrix  using  the  Cholesky  method.  The 
manipulations  are  carried  out  by  an  especially  efficient  routine 
developed  by  M.  Vander  Vorst  of  JAYCOR. 
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1 


stability  of  the  solution  requires  that  the  time  step  does  not 
exceed  values  imposed  by  the  advection  and  diffusion  terms. 


(St  < min 


6x 


fiy. 


1.  J 


i.J 


R ^ X \ 2 

4 («x) 


f 


2 


(31) 


Experience  has  shown  that  if  6t  is  taken  to  be  one-quarter  of 
the  maximum  time  step  allowed,  then  at  least  three  significant 
figures  are  obtained  in  the  velocity  and  pressure  fields  with 
five  iterations. 


4.  COMPARISON  WITH  LINEAR  STABILITY  THEORY 


The  linearized  form  of  equations  (13)  - (15)  have  been  studied 
in  detail,  e.g.,  Wazzan,  et  al . , (1968),  Jordinson  (1970),  the 
case  of  greatest  interest  being  when  U takes  on  the  Blasius 
profile.  See  also  the  reviews  by  Tani  (1969)  and  Reshotko 
(1976).  Not  only  have  detailed  quantitative  solutions  been 
obtained,  but  the  predicted  fluid  dynamic  behavior  is  in  close 
agreement  with  experimental  observation  of  small  disturbances 
in  the  boundary  layer  of  a flat  plate.  For  example,  see  Ross, 
et  al.,  (1970).  Comparison  with  LST  offers  then  a meaningful 
test  of  the  numerical  approach. 

When  the  disturbances  are  periodic  in  space,  the  eigenfunctions 
are  characterized  by  the  two  parameters  oc  and  R.  Figure  3 shows 
the  curve  of  neutral  stability,  that  is  those  distrubances  which 
neither  grow  nor  die  out,  as  determined  by  Jordinson;  points 
inside  the  curve  are  unstable  flows,  those  outside  are  stable. 
Below  the  critical  Reynolds  number  of  520,  LST  predicts  that  all 
disturbances  die  out  and  thus  that  the  boundary  layer  is 
completely  stable. 

To  compare  with  the  predictions  of  LST  a series  of  numerical 
experiments  were  conducted.  In  each  case  the  mesh  consisted  of 
18  constant-sized  cells  in  the  horizontal,  periodic  direction 
and  34  variably-sized  cells  in  the  vertical.  The  smallest 
vertical  cells  were  0.15  wide,  while  above  y = 2.06  the  cells 
expanded  geometrically  so  that  the  top  of  the  last  cell  was 
almost  1005  from  the  solid  boundary.  Figure  1 shows  the  mesh 
configuration  up  to  y = 255.  The  top  boundary  appears  to  not 
influence  the  results,  since  placing  it  further  away  led  to 
no  change  in  the  flow  properties.  The  base  flow  was  given  a 
small  disturbance  which  was  then  allowed  to  evolve  toward  a 
final  state. 
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Figure  3.  Neutral  Stability  Curve  and  Cases  for  which  Numerical  Calculations  have 
been  made.  Numbers  indicate  the  logarithmic  energy  change  per  period. 


The  circles  shown  in  Figure  3 represent  the  matrix  of  numerical 
experiments  conducted.  The  numbers  indicate  the  logarithmic 
energy  change  per  period  of  oscillation 


0 = Jin  ^ > T = period  of  disturbance.  (32) 

This  quantity  has  not  yet  reached  a final  value  in  all  cases, 
so  the  quoted  value,  taken  at  the  end  of  a run,  may  change 
slightly.  Figure  4 shows  the  variation  with  time  of  the  energy 
for  most  of  the  runs. 

Figure  5 compares  the  time  variation  of  production  and  dissi- 
pation for  one  case  near  the  neutral  curve,  where 

production  = JT  dx  dy  ( - uv  (33) 

2 2 

dissipation  = -JXdx  dy  (uV  u + vV  v)  . (34) 

The  production  is  the  rate  at  which  energy  is  extracted  from 
the  base  flow,  while  dissipation  is  the  rate  at  which  it  is  lost 
to  heat.  The  rates  should  be  equal  for  a neutral  case.  Since 
the  initial  disturbance  is  not  a pure  eigenmode  of  the  system, 
and  since  neighboring  modes  die  out  slowly,  the  flow  requires 
a considerable  time  to  settle  into  a single,  pure  motion.  That 
difference  distinguishes  a true  initial  value  problem  from  the 
eigenvalue  problem  posed  by  LST. 

Since  negative  values  in  Figure  3 indicate  damped  motion  and 
positive  values  indicate  growth,  the  results  are  not  in  conflict 
with  LST.  In  fact,  when  nearby  pairs  of  numerical  calculations, 
one  showing  growth  and  the  other  decay,  are  used  to  interpolate 
a neutral  condition,  the  results  agree  with  LST  to  within  a few 
percent  for  both  the  wavenumber  and  frequency,  see  Table  1. 

Here,  0 is  the  dimensionless  frequency  of  oscillation.  The 
numerical  calculation  appears  to  have  correctly  reproduced  the 
overall  energy  balance  and  time  behavior  in  the  flow. 
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Figure  5.  Comparison  of  Production  and  Dissipation 
Rates  During  the  Relaxation  to  a Final 
State  for  a Case  Near  the  Neutral  Curve. 
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Table  1. 


1000 


2000 


2600 


3000 


a 

B 

numerical 

0.348 

0.128 

LST 

0.351 

0.131 

numerical 

0.1239 

0.0327 

LST 

0.1234 

0.0332 

numerical 

0.295 

0.0927 

LST 

0.297 

0.0947 

numerical 

0.1051 

0.0257 

LST 

0.1029 

0.0235 

0.30 


numerical 

LST 


0.1167 
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Next,  details  of  the  flow  field  were  examined.  Figure  6 
shows  the  variation  with  time  of  the  horizontal  disturbance 
velocity  at  a particular  streamwise  position  and  for  various 
distances  from  the  wall.  This  format  is  analogous  to  the 
continuous  output  of  a velocity  probe  and  shows  the  experi- 
mentally observed  pattern  (Schlichting,  1969,  p.403)  that  the 
amplitude  grows  with  distance  from  the  wall,  reaches  a peak, 
suffers  a phase  reversal,  and  dies  out.  The  calculated  period 
of  oscillation  is  93.79  while  the  theoretical  value  is  93.48. 

Figures  (7)  and  (8)  compare  the  amplitude  and  phase  variation 
with  distance  from  the  wall  found  from  a numerical  calculation 
and  Jordinson's  eigenfunction  solution.  The  comparison  shows 
that  the  details  of  the  flow  are  reproduced  completely, 
including  the  subtle  phase  variations  occurring  in  the  flow  near 
the  wall.  The  numerical  calculation  appears  to  have  correctly 
reproduced  the  flow  details  as  well  as  the  overall  energy 
balances . 


Figure  6.  Time  Traces  of  the  Horizontal  Velocity 
Component  for  Various  Distances  from 
the  Wall. 
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Figure  7.  Comparison  of  Horizontal  Velocity 
Amplitude  Variation. 
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5.  FINITE  AMPLITUDE  EFFECTS 


The  previous  comparisons  with  linear  stability  theory  suggest 
that  the  behavior  of  small  amplitude  disturbances  is  being 
correctly  computed  numerically.  A series  of  calculations  was 
then  performed  to  determine  the  effects  of  finite  amplitude. 

In  each  case  the  initial  disturbance  had  the  same  spatial  dis- 
tribution but  the  maximum  amplitude  of  u'  was  taken  as  0.0001, 
0.01,  0.03,  and  0.1  of  U^.  The  results  show  a trend  conjectured 
by  nonlinear  stability  theories  (Stewartson,  1975)  that  finite 
amplitude  effects  can  be  destabilizing. 

Figure  9 shows  the  variation  of  the  total  energy  in  the 
disturbance  flow  in  each  case,  normalized  to  the  same  initial 
value.  The  case  considered  is  R = 500,  a = 0.30,  a point  in 
the  subcritical  regime  where  LST  predicts  all  disturbances  to 
die  out.  At  the  smallest  amplitudes  the  disturbance  dies  out 
at  a rate  consistant  with  LST.  For  u'  = 0.03,  however,  the 
decay  has  been  noticeably  decreased,  while  for  u'  = 0.1  the 
disturbance  is  almost  neutrally  stable.  Although  a truly 
unstable  state  was  not  achieved,  the  destabilizing  effects 
seem  to  be  definite.  A similar  variation  is  seen  for  unstable 
small  amplitude  modes:  as  they  grow  to  finite  amplitude  the 
rate  of  growth  increases. 

Figure  10  further  illustrates  the  effects  of  finite  amplitude, 
by  considering  the  integrated  production  and  dissipation  rates 
for  the  case  R = 500,  a = 0.30.  At  small  amplitudes  dissipation 
always  exceeds  production  and  the  disturbance  monotonically  dies 
out.  At  u'  = 0.1,  however,  the  production  exceeds  dissipation 
during  part  of  the  oscillation,  so  that  for  some  time  it  grows 
and  for  some  time  it  is  damped.  The  damped  periods  are  greater 
so  that  the  disturbance  is  dying  out  overall,  but  a balance  is 
almost  being  struck. 


The  accuracy  of  finite  difference  techniques  for  large 
amplitude  phenomena  is  difficult  to  assess  and  so  a detailed 
investigation  is  probably  not  warranted  at  this  time.  Still, 
the  results,  indicate  that  the  numerical  solution  technique 
used  in  this  work  may  be  a useful  tool  in  studying  finite 
amplitude  effects  in  the  boundary  layer. 
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6 . CONCLUDING  REMARKS 


The  finite-difference  method  presented  here  has  been  shown  to 
accurately  reproduce  analytical  results  concerning  boundary 
layer  stability  with  a modest  storage  requirement  (612  interior 
cells)  and  modest  computational  time  (typically  2-5  minutes 
of  CDC  7600  central  processor  time).  This  efficiency,  coupled 
with  encouraging  results  in  the  nonlinear  regime,  makes  the 
approach  an  attractive  candidate  for  studying  the  large 
amplitude,  three-dimensional  motion  that  accompanies  boundary- 
layer  transition. 
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